R Notebook

Install/Load Required Packages

We’ll use the tidyverse. In addition, we’ll use arule to implement association rules mining.

packages <- c("tidyverse", "haven", "arules", "plotly")

# Installing Packages
# install.packages(packages, update = TRUE, ask = FALSE) 
  
# Loading Packages
for(package in packages) {
  do.call("library", list(package))
}

Data Processing

We’ll first load the data from the NHANES website (data documentation) to a tibble. We’ll only look at prescriptions taken within the last month, and we’ll only need the patient identifier and drug name columns.

# Load data from NHANES into prescriptions_df
prescriptions_df <- read_xpt("https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/RXQ_RX_J.XPT")

# Filter to prescriptions that were taken in the last month
# Select patients (SEQN) and drug name (RXXDRUG) 
prescriptions_df <- prescriptions_df %>%
  filter(RXDUSE == 1) %>% 
  select(SEQN, RXDDRUG)

# View dataframe
prescriptions_df

For arules to work, we’ll need to group the data by the patient identifier, and we’ll the drug names for each patients to be stored within a list in the tibble. As a final step, we removed all patients with no drugs.

prescriptions_df <- prescriptions_df %>% 
  filter(
    RXDDRUG != "55555",
    RXDDRUG != "77777",
    RXDDRUG != "99999",
    RXDDRUG != ""
  ) %>%
  group_by(SEQN) %>% 
  summarize(RXDDRUG = list(unique(RXDDRUG)))

prescriptions_df

Mine Frequent Itemsets

Now we can implement association rules mining. First we’ll mine the frequent itemsets using the apriori algorithm. We’ll generate a list of itemsets with a support of at least 0.005.

frequent_itemsets <- apriori(
  as(prescriptions_df$RXDDRUG, "transactions"), 
  parameter = list(supp=0.005, target="frequent itemsets")
)

Apriori

Parameter specification:

Algorithmic control:

Absolute minimum support count: 19 

set item appearances ...[0 item(s)] done [0.00s].
set transactions ...[666 item(s), 3951 transaction(s)] done [0.00s].
sorting and recoding items ... [140 item(s)] done [0.00s].
creating transaction tree ... done [0.00s].
checking subsets of size 1 2 3 4 done [0.00s].
sorting transactions ... done [0.00s].
writing ... [360 set(s)] done [0.00s].
creating S4 object  ... done [0.00s].

frequent_itemsets <- inspect(frequent_itemsets) %>% 
  arrange(desc(support))

frequent_itemsets

Next let’s look at association rules. From looking at the rules with the highest lift, we see that the insulin drugs are highly associated. What are other patterns that may be interesting?

rules <- apriori(
  as(prescriptions_df$RXDDRUG, "transactions"), 
  parameter = list(supp=0.005, conf=0.005, target="rules")
)

Apriori

Parameter specification:

Algorithmic control:

Absolute minimum support count: 19 

set item appearances ...[0 item(s)] done [0.00s].
set transactions ...[666 item(s), 3951 transaction(s)] done [0.00s].
sorting and recoding items ... [140 item(s)] done [0.00s].
creating transaction tree ... done [0.00s].
checking subsets of size 1 2 3 4 done [0.00s].
writing ... [618 rule(s)] done [0.00s].
creating S4 object  ... done [0.00s].

rules <-data.frame(inspect(rules)) %>%
  filter(
    lift >= 1,
    lhs != "{}",
    rhs != "{}"
  ) %>%
  select(-Var.2) %>%
  arrange(desc(lift))

rules

rules %>%
  filter(
    support > 0.02 &
    confidence > 0.6
  ) %>%
  select(
    lhs,
    rhs
  )

NA

Visualize Rules

We can plot the support and confidence (or lift) or each rule on a scatterplot.

ggplot(rules) +
  aes(x=support, y=confidence, label=lhs, label2=rhs) +
  geom_point()

ggplotly()